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We performed calculations of electronic, optical and transport properties of graphene on hBN with 
realistic moire patterns. The latter are produced by structural relaxation using a fully atomistic 
model. This relaxation turns out to be crucially important for electronic properties. We describe 
experimentally observed features such as additional Dirac points and the ’’Hofstadter butterfly” 
structure of energy levels in a magnetic field. We find that the electronic structure is sensitive to 
many-body renormalization of the local energy gap. 
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The physical properties of van der Waals heterostruc¬ 
tures can change drastically in comparison with the ones 
of the constituent two-dimensional materials pQ. Re¬ 
cent experiments of graphene on hexagonal boron-nitride 
(hBN) show that hBN can act like an effective peri¬ 
odic potential for graphene, leading to secondary Dirac 
points 0 0. The graphene/hBN heterostructures are 
of fundamental interest as an example of a quantum me¬ 
chanical system with tunable incommensurate potentials. 
Such incommensurate potentials are important for qua¬ 
sicrystals [4 but they are not tunable, whereas in the 
graphene-hBN systems it is possible to change the po¬ 
tential by changing the mutual orientation of graphene 
and hBN layers. It was long predicted that a system 
under influence of both a crystal potential and a mag¬ 
netic field, with a magnetic period incommensurate with 
that of the crystal, would exhibit a recursive spectrum 
now called Hofstadter’s butterfly [5], which has been ob¬ 
served in experiments with misaligned graphene on hBN 
in 2013 [SHE]- 

Although hBN has a structure similar to that of 
graphene, the lattice mismatch of 1.8% will cause moire 
patterns, meaning that there is no uniform stacking in 
the sample. An extra difficulty is posed by the recently 
observed transition at very small angles from an incom¬ 
mensurate state, with little deformation of graphene, to a 
commensurate state, where regions of stretched graphene 
are separated by domain walls 0. The computational 
challenge lies in the fact that at such angles these su¬ 
perlattices have unit cells consisting of tens of thousands 
atoms, making it impossible to study them using meth¬ 
ods such as DFT. The tight-binding propagation method 
(TBPM) described in [TO] . can be used to study systems 
with hundreds of millions of atoms, circumventing this 
problem. In this research, we apply TBPM to graphene- 
hBN heterostructures with various rotation angles, based 
on structures determined by atomistic simulations of re¬ 
alistic moire patterns. 


Various approaches have been developed to describe 
graphene-hBN with an effective tight-binding (TB) 
Hamiltonian mm- While these methods give reason¬ 
able results, they lack the flexibility needed to apply them 
to other systems. We present an approach consisting of 
three parts, namely: i) structural relaxation of graphene 
on top of hBN with an empirical potential, ii) modifi¬ 
cation of the TB parameters due to this relaxation, and 
iii) calculation of electronic properties with these mod¬ 
ified TB-parameters. There are multiple advantages to 
this approach. First, the construction of the TB model is 
solely based on the three-dimensional coordinates of the 
carbon atoms, thus one could use this same method for 
graphene on top of other substrates or graphene under 
mechanical strain. Second, it is easy to incorporate extra 
disorder such as carbon vacancies, ad-atoms and ripples, 
etc. 

The first step is the relaxation of graphene on hBN. 
We follow the approach of Ref. m , where it was shown 
that moire patterns can be used as a probe of interpla- 
nar interactions for graphene on hBN. We construct su¬ 
percells of rotated graphene on hBN with misorientation 
angles 6 and corresponding moire patterns with period 
A 0 . The graphene atoms interact through the reac¬ 
tive empirical bond order potential REBO m, as imple¬ 
mented in the molecular dynamics code LAMMPS 16 . 
The hBN substrate is kept rigid, mimicking a bulk sub¬ 
strate. As no empirical potential for the interactions be¬ 
tween graphene and hBN is available, we use the registry- 
dependent Kolmogorov-Crespi potential m developed 
for graphite. We neglect the correction for bending in¬ 
troduced to describe carbon nanotubes. We set the ra¬ 
tio of C-B/C-N interactions to 30% with the C-N inter¬ 
action twice as strong as the original C-C interaction, 
as this leads to better agreement with experimental re¬ 
sults and ab initio calculations 00. We mini¬ 

mize the total potential energy by relaxing the graphene 
layer by means of FIRE [20] . a damped dynamics algo- 
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FIG. 1. The modified TB parameters for a relaxed sample of graphene on hBN with 0 — 0° (A = 13.8 nm). From left to right 
the on-site potential v, and the hopping parameters G, t <2 and 1 3 . The color bars are in units of t — 2.7 eV. 


rithm. For aligned samples (0 = 0°), this relaxation leads 
to significant changes in bond length along the moire pat¬ 
tern. The degree of deformation decreases with increas¬ 
ing angle [2T] . 

After relaxation, we use the following graphene TB- 
Hamiltonian. The main idea of our method is that the 
TB parameters are modified as a function of a small dis¬ 
placement out of equilibrium of the carbon atoms. The 
general TB-Hamiltonian for graphene is given by: 

H = -J2 Vic\cj, ( 1 ) 

(i,j) i 

where only the nearest-neighbor hopping and on-site po¬ 
tential are taken into account. Including next-nearest- 
neighbor hoppings will result in minor changes [22]. The 
change in the hopping parameter tij can be written 
as [23] : 



(a) (b) 

FIG. 2 . (a) Numerical results for the DOS of unrelaxed and 
relaxed graphene, (b) DOS for different angles 0. The ex¬ 
tra cones move outward, indicated by the arrows, and disap¬ 
pear for large angles. The corresponding moire lengths A are: 
13.8 nm, 11.9 nm and 6.7 nm respectively. 


Uj = t exp(—3.37(r,j/ao - 1)), (2) 

where is the distance between atoms i and j, t = 
2.7 eV is the regular hopping parameter, and ao = 
1.42 A is the equilibrium carbon-carbon distance for 
graphene. For the on-site potential Vi we calculate an 
effective area S of each carbon atom [24], which will be 
changed due to local deformations resulting in a modu¬ 
lated value for vf. 


where gi = 4 eV. This value corresponds to the screened 
deformation potential, which gives reasonable description 
of transport properties m, and is close to density func¬ 
tional estimates [26]. Figure [l] shows the change of the 
on-site potential and of the hopping parameters for a 
relaxed layer of graphene on hBN with rotation angle 
0 = 0°. A clear periodic modulation with period A is 
found in all parameters. 

Electronic properties are calculated using the TBPM, 
a method based on the propagation of the wavefunc- 
tion of the time-dependent Schrodinger equation using 
Chebychev polynomials mm- The correlation func¬ 
tion (</>(0)|e -zi ^|</>(0)} is calculated at each time step. 
The density of states (DOS) can then be obtained by a 


Fourier transform of these correlation function. To in¬ 
crease the accuracy of the electronic calculations the su¬ 
percells are repeated so that the total system consists of 
~ 6000 x 6000 carbon atoms. 

The first step to validate our method is to compare the 
DOS of pristine graphene (unrelaxed) to that of graphene 
on hBN after energy minimization (relaxed), as shown in 
Figure [2a] Secondary Dirac cones appear at both the 
electron and hole side, as seen in experiments |2 El [SHE]. 
The positions depend on the reciprocal lattice vector G of 
the moire pattern, and is given by Ed = ±hvp |G| /2 = 
-E2'kHvf /(a/3A) where vp is the Fermi velocity [llj 28, 
29[. Due to the substrate the depth of the extra cones is 
asymmetric and highly dependent on the value of the on¬ 
site potential. The position of the extra Dirac cones will 
change with misorientation angle 6 as A depends on 9M- 
Figure [2b] shows how small angular variations shift the 
extra cones. The effect of the relaxation decreases with 
increasing #, meaning that the differences of the DOS 
also become negligible for large 6. 

The real-space distribution of eigenstates can be com¬ 
pared with the LDOS images obtained from STM mea¬ 
surements. In general it is hard to obtain the eigen¬ 
states corresponding to a TB Hamiltonian of a system 
with millions of atoms. We obtain the so called quasi- 
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FIG. 3. Amplitude of the quasi-eigenstates for different ener¬ 
gies for 6 — 0°. The left-hand panels show sublattice A and 
the right-hand panels show sublattice B. For energies closer 
to the extra Dirac cones a clear moire pattern can be distin¬ 
guished. Only roughly one thousands of the system is shown. 


eigenstates m , which are close to the real eigenstates, 
by using the TBPM. Figure [3] shows the amplitude of 
some of these quasi-eigenstates close to the additional 
Dirac cones. The hBN substrate breaks the sublattice 
symmetry, and therefore we plot the quasi-eigenstates 
separately. Some localization is found for the quasi¬ 
eigenstates. We see that for energies close to the Fermi 
energy the difference between amplitudes is negligible. 
For energies closer to the additional Dirac cones a clear 
moire pattern can be distinguished. 
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FIG. 4. Density of states (left) and DC conductivity a (right) 
as a function of the carrier density n e for 6 — 0° and for 
varying A U. 

The appearance of additional Dirac cones in the DOS 
and the signatures of localization in the quasi-eigenstates 
indicate that the electronic structure is strongly influ¬ 
enced by relaxation. The transport measurements of 


DC conductivity of graphene on hBN in recent experi¬ 
ments [SHE] show clearly asymmetric drops of the con¬ 
ductivity at the secondary Dirac points on the hole and 
electron sides. The decreasing of the conductivity on 
the hole side is more significant, with a value even lower 
than the minimum conductivity at the Dirac point in 
Ref. [7. We calculate the DC conductivity by using the 
Kubo formalism [30 within the TBPM [10]. The results 
shown in figure [4] do not have such minimum on the hole 
side as in experiments. It could be obtained by using 
interaction strength in the empirical potential used for 
the relaxation m much stronger than is suggested by 
ab initio total energy calculations m- However, there is 
an interaction which we have not yet considered, namely 
the local gap opening induced by the substrate [IS]. It 
is known that the many-body effects can increase the 
gap dramatically [32 , and more accurate GW calcula¬ 
tion gives a several times larger gap m in comparison 
with DFT [18]. To take into account the sublattice asym¬ 
metry due to many-body effect, we add a local gap term 
according to the potential difference between one site and 
its three neighbors as: 


/ . A .92 

Vi = Vi + A Vi = Vi + — 


o F V i+S 


8 = 1,2,3 


(4) 


The strength of the local gap, which is controlled by 
the parameter g 2 in Eq. [4] is given by the average of 
the potential difference between sublattices A and B 
A U = (|A^|). Numerical calculations of the DOS in 
Figure [4] show that the depth of the additional minima 
at energy Ejj can be tuned by the local gap AC/. For in¬ 
creasing AC/, the minimum on the hole side of the DOS 
becomes deeper, while the one on the electron side first 
disappears for small AU and then reappears for large 
AC/. Although it is very difficult to estimate A U accu¬ 
rately since there is no quantitatively accurate theory of 
many-body effects in graphene, we can use the one ob¬ 
tained by Bokdam et al. [19], a GW band gap of 32 meV 
for incommensurable graphene on hBN with 0 = 0° as a 
reference value. For A U = 32 meV, we see clearly a de¬ 
crease (increase) of DOS and DC conductivity (a) at the 
extra Dirac point on the hole (electron) side. The trans¬ 
port calculation with AU = 32 meV reproduces well the 
experimental observations in Ref. m n]. On the other 
hand, in Ref. [7], the value of DC conductivity at the 
extra Dirac point on the hole side is smaller than the 
minimum conductivity at the Dirac point. This is only 
possible by using a larger AC/, for example, a drops to 
zero by doubling AU as 64 meV. Our numerical results 
suggest that the experimentally observed insulating state 
at the extra Dirac point on the hole side |HE] is a signa¬ 
ture of strong local gap induced by many-body effects. 

In the presence of a perpendicular magnetic field, 
the quantization of the energy eigenstates leads to dis- 
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crete Landau levels. The modulation induced by the 
moire patterns, splits the flat Landau bands of pristine 
graphene into minibands, the so called ’’Hofstadter but¬ 
terfly spectrum” which has been conformed in several re¬ 
cent experiments [6H8| . In order to verify the splitting of 
the Landau levels in our TB model, we show the contour 
plot of DOS as a function of magnetic field strengths in 
Figure |5j For both A U = 32 and A U = 64 meV, there is 
a clear splitting of the Landau levels with increasing mag¬ 
netic field, and the splitting becomes more clear when the 
stronger local gap term is included. 




E/t 

FIG. 5. The density of states for varying magnetic field with 
(top) A U — 32 meV and (bottom) A U = 64 meV. Strong 
extra peaks at both electron and hole side are observed, which 
for higher magnetic fields split into two. 




FIG. 6. Optical conductivity a for various angles 6 with (left) 
A U = 32 meV and (right) A U = 64 meV. ao = ire 2 /2h is the 
universal optical conductivity of graphene. Inset: the optical 
conductivity for larger scales. 


Another quantity of great experimental and practical 
interest is the optical conductivity, that we calculate by 
using the TBPM [lTJ, 33 . Due to the presence of moire 
pattern, we expect that there should be signatures of the 
extra Dirac cones in the optical spectrum. Figure [6] shows 
the optical spectrum a of graphene on hBN with three 
different orientation angles 6. For high energies the en¬ 
hanced peak around uj = 2£, resulting from the optical 
transition between Van Hove singularities at E = d=t, is 
similar to pristine graphene. Futhermore, there are addi¬ 
tional peaks at photon energy about u) = 2 \E&\ (around 
0.1 ^ 0.2 t, depending on the angle 0 ), correspond¬ 
ing to the optical transitions between the peak states 
around the extra Dirac points on the hole and electron 
sides. The amplitudes of these peaks increase signifi¬ 
cantly with larger local gap term. It is known that the 
optical conductivity of graphene for visible light has a 
universal value, our results with moire patterns indicate 
that the optical conductivity becomes tunable by chang¬ 
ing the relative orientations between graphene and its 
hBN substrate. 

To conclude, we have shown that merely taking into ac¬ 
count the periodic modulation in graphene caused by a 
substrate is enough to describe new features in the elec¬ 
tronic and optical properties of graphene. The many- 
body enhancement of the local energy gap is crucially 
important to reproduce the experimentally observed in¬ 
sulating state at the extra Dirac point on the hole side. 
We also show that the optical conductivity of graphene 
is tunable by varying the relative orientations between 
graphene and its hBN substrate. The presented ap¬ 
proach for the construction of the TB model is not lim¬ 
ited to graphene-BN heterostructures, but can be used 
for graphene with other substrates, such as Ru and Cu, 
and can be extended to include various types of disorder. 
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